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ABSTRACT 

An  acoustic  tomography  array  consisting  of  six  transceiver  moorings  was 
jointly  deployed  by  Woods  Hole  Oceanographic  Institution  and  Scripps 
Institution  of  Oceanography  in  the  Greenland  Sea  during  the  second  half  of  1988. 
Two  of  the  primary  objectives  of  this  thesis  are:  (1)  to  set  up  and  test  a  stochastic 
3-D  inversion  code  for  the  Greenland  Sea  Acoustic  Tomography  data  analysis; 
and  (2)  to  evaluate  the  performance  of  the  acoustic  system  through  resolution 
and  variance  analyses.  In  acoustic  tomography,  the  sound  speed  perturbation 
field  is  estimated  from  measured  acoustic  travel  time  perturbation  data.  A  unique 
sound  speed  perturbation  estimate  can  be  constructed  using  the  Guass-Markoff 
theorem.  However,  the  theorem  requires  the  specification  of  the  covariance  of 
the  sound  speed  perturbation  field,  which  is  generally  not  exactly  known.  Via 
computer  simulation,  we  examined  the  sensitivity  of  the  estimate  to  uncertainty 
in  the  sound  speed  field  correlation  specified.  In  addition,  we  also  examined  the 
effects  of  an  increased  random  experimental  noise  level  and  a  change  in  array 
geometry  due  to  mooring  failure  on  the  estimate.  The  three  major  results  are 
that:  (1)  the  estimate  is  less  sensitive  to  a  positive  uncertainty  in  correlation 
length  than  to  a  negative  uncertainty  in  an  ocean  volume  containing  large 
structures,  while  it  is  more  sensitive  to  a  positive  uncertainty  than  to  a  negative 
uncertainty  in  an  ocean  volume  containing  small  structures;  (2)  the  estimate 
error  is  primarily  bias  error  rather  than  random  error;  and  (3)  the  failure  of  a 
mooring  causes  a  large  increase  in  RMS  error  in  regions  no  longer  containing 
acoustic  rays,  but  it  results  in  an  increase  in  RMS  error  of  only  25%  in  regions 
which  still  contain  acoustic  rays. 
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I.    INTRODUCTION 

A.   OCEAN  ACOUSTIC  TOMOGRAPHY 

Ocean  Acoustic  Tomography  is  a  method  used  to  monitor  the  mesoscale 
ocean  variability  (which  is  the  oceanic  analog  of  atmospheric  weather)  and  it  was 
introduced  by  Munk  and  Wunsch  (1979).  This  technique  is  analogous  to  the 
medical  X-ray  procedure  known  as  Computer  Assisted  Tomography  (CAT) 
(Figure  1-1. a).  Roughly  speaking,  tomography  exploits  the  fact  that  the  ocean  is 
"transparent"  to  acoustic  rays  to  remotely  sense  the  properties  of  an  ocean 
region. 


body 


X-ray  source 


ocean  front 


acoustic 
transceiver 


(a)  Medical  CAT 


(b)  Ocean  Acoustic  Tomography 


Figure  1-1:  The  Comparison  of  Medical  CAT  and  Ocean  Acoustic 
Tomography. 

In  practice,  a  number  of  acoustic  transceivers  are  deployed  at  positions  chosen  to 

allow  for  coverage  of  an  ocean  volume  of  interest  (such  as  a  region  containing 

mesoscale  eddies  or  a  frontal  system)  (Figure  1-l.b).  The  most  common 


application  of  tomography  is  for  estimating  the  perturbation  of  the  sound  speed 
field  from  a  set  of  measured  acoustic  travel  time  perturbations.  The 
perturbations  in  sound  speed  are  assumed  to  be  so  small  that  the  perturbations  in 
acoustic  travel  time  between  each  pair  of  transceivers  are  linearly  related  to  the 
sound  speed  perturbations.  The  modeling  of  the  travel  time  perturbations  due  to 
the  sound  speed  perturbations  is  known  as  the  forward  problem.  Once  the 
forward  problem  is  solved,  inverse  methods  which  are  widely  used  in 
geophysical  research  (Backus  and  Gilbert,  1967)  are  applied  to  the  travel  time 
data  for  the  reconstruction  of  the  the  sound  speed  perturbation  field. 

Ocean  Acoustic  Tomography  offers  several  advantages  over  conventional 
hydrographic  surveying  method.  These  advantages  are  pointed  out  by  Chiu 
(1978):  (1)  the  system  can  be  implanted  in  the  ocean  on  a  semipermanent  basis  to 
allow  for  continuous  observation;  (2)  it  is  not  affected  greatly  by  weather 
conditions;  (3)  it  has  high  temporal  resolution;  (4)  it  can  cover  an  extensive 
volume  of  the  ocean  interior  and  probe  the  different  parts  simultaneously;  and 
(5)  only  a  few  moorings  are  needed,  thus  minimizing  the  effort  in  deployment 
and  maintenance. 

Since  the  first  successful  experiment  (the  1981  Three-dimensional  Mesoscale 
Experiment),  additional  tomography  projects  have  provided  measurements  of 
mesoscale  eddies  (Cornuelle  et  al,  1983),  planetary  waves  (Chiu  et  al,  1987), 
currents  (DeFerrari  et  al,  1986),  internal  waves  (Stoughton  et  al,  1986),  basin 
mode  oscillations  (Bushong,  1987),  and  surface  waves  (Lynch  et  al,  1987).  In  the 
future,  monitoring  of  large-scale  ocean  dynamics  on  a  global  basis  may  be 
achieved  using  cross-basin  transmissions. 


B.    GREENLAND  SEA  PROJECT  OCEAN  ACOUSTIC 
TOMOGRAPHY 

The  Greenland  Sea  Project  (GSP)  is  a  plan  developed  by  the  international 
Greenland  Sea  Science  Planning  Group  which  was  appointed  by  the  Arctic  Ocean 
Sciences  Board  (AOSB).  The  overall  goal  of  this  five-year  program  (from  1987 
to  1992)  defined  by  AOSB  is  to  understand  the  large  scale,  long-term 
interactions  among  the  air,  sea,  and  ice  in  the  Greenland  Sea.  The  primary 
region  of  the  study  is  bounded  by  Fram  Strait  to  the  north,  Spitsbergen  and  the 
Mohn  Rise  to  the  east,  the  Greenland-Jan  Mayen  Ridge  to  the  south,  and 
Greenland  to  the  west  (Greenland  Sea  Science  Planning  Group,  1986,  pp.  1-7). 

The  plan  is  designed  to  study  the  following  ocean  dynamics:  (1)  the  seasonal 
and  interannual  variability  of  the  sea  ice  cover;  (2)  ocean  ventilation  and 
convection  of  the  deep  water;  (3)  ocean  circulation  and  mixing;  (4)  atmosphere 
energetics;  and  (5)  biological  processes.  The  Ocean  Acoustic  Tomography  Array 
is  a  component  used  in  GSP  to  monitor  the  process  of  ocean  ventilation  and 
convection  in  the  Greenland  sea  central  gyre  (Greenland  Sea  Science  Planning 
Group,  1986,  pp.  1-7).  It  is  the  process  of  ventilation  and  convection  that  gives 
the  Greenland  Sea  central  gyre  the  ability  to  affect  the  oceans  throughout  the 
world. 

An  acoustic  tomography  array  consisting  of  six  transceiver  moorings  with  a 
pentagonal  geometry  was  jointly  deployed  by  Woods  Hole  Oceanographic 
Institution  and  Scripps  Institution  of  Oceanography  in  the  Greenland  Sea  during 
the  second  half  of  1988.  Figure  1-2  shows  the  tomography  array  position  and  the 
GSP  area. 
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Figure  1-2:  Acoustic  Tomography  Mooring  Array  Position  and  Geometry 
Configuration. 


For  convenience  of  calculation,  the  given  geodetic  position  (latitude  and 
longitude)  of  each  array  element  has  been  translated  into  an  xy  position  using  a 
Transverse  Mercator  (TM)  projection.  To  do  this,  the  longitude  of  the  central 
mooring  (array  #6),  was  chosen  as  the  central  meridian  of  the  TM  projection. 
The  false  origin  was  set  at  73°N  and  shifted  to  120  km  west  of  the  central 
meridian.  This  coordinate  system  gives  only  slight  position  distortion  at  the  edge 
of  the  array.  The  ray  paths  can  then  be  calculated  in  the  xy  planar  coordinates 
rather  than  in  geodetic  coordinates.  TABLE  1-1  shows  the  coordinate  conversion 
for  the  moorings  as  well  as  the  depths  of  the  acoustic  sources  and  vertical 
receiver  arrays. 

TABLE  1-1:  COORDINATE  CONVERSION  OF  THE  ACOUSTIC 
MOORINGS  (USING  WGS72  SPHEROID) 


Mooring         Lat. 


Lon. 


1  75°47.9'N  1°04.7,W 

2  74°53.3'N  1°24.8'E 

3  74°00.4'N  1'05.2'W 

4  74°18.9'N  4'59.0'W 

5  75'25.2'N  5°15.9*W 

6  74°54.0*N  2°12.0'W 


X 

(meter) 

y 

(meter) 

S-depth* 
(meter) 

R-depth 
(meter) 

150726 

312598 

94.5 

145.3 

225124 

213940 

95.0 

145.4 

154249 

112659 

94.6 

145.1 

36008 

148715 

94.6 

117.0 

33882 

272309 

94.8 

117.7 

120000 

212040 

95.2 

145.7 

*   The  depth  of  acoustic  source. 

**  The  depth  of  vertical  acoustic  receiver  array. 


C.   THESIS  OBJECTIVE 

The  purpose  of  this  research  is  to  develop  computational  tools  in  preparation 
for  analysis  of  the  tomography  data  from  the  GSP.  Specifically,  our  objectives 
are: 

•  to  set  up  and  test  a  stochastic  3-D  inversion  code  for  the  Greenland  Sea 
Acoustic  Tomography  data  analysis;  and 

•  to  investigate  the  sensitivity  of  the  sound  speed  perturbation  estimate  to 
our  uncertainty  in  the  sound  speed  field  correlation,  changes  in  the  random 
experimental  noise  level,  and  changes  in  array  geometry  due  to  mooring  failure. 
The  second  objective  is  accomplished  via  computer  simulation. 

In  Chapter  II,  we  discuss  the  behavior  of  an  acoustic  ray  in  an 
inhomogeneous  moving  medium.  A  linear  sound  speed  profile  is  taken  as  the 
reference  state  of  the  sound  speed  field.  Based  on  this  reference  state,  a 
numerical  4th  order  Runge-Kutta  integration  method  was  used  to  calculate  the 
paths  of  the  eigenrays  for  establishing  the  forward  acoustic  model.  We  then 
discretize  the  forward  model  by  dividing  the  ocean  volume  into  500  boxes  in 
order  to  cast  the  problem  into  a  matrix  form  for  the  computer  simulation. 

In  Chapter  III  we  present  a  three-dimensional  stochastic  inverse  method 
which  is  the  distribution-free  Gauss-Markoff  estimator.  Due  to  insufficient 
experiment  data,  the  inverse  problem  is  underdetermined.  In  this  stochastic 
approach,  a  priori  information  is  specified  in  the  covariance  matrix  of  sound 
speed  perturbations.  The  covariance  gives  additional  constraints  to  the  system 
and  therefore  a  unique  solution  is  obtained.  Two  measures,  RMS  error  and 
resolution  length,  are  used  to  quantify  the  performance  of  the  estimator  at  each 
box  location. 


The  ocean  is  a  dynamic  and  inhomogeneous  environment.  It  is  difficult  to 
incorporate  an  exact  sound  speed  perturbation  covariance  matrix  as  a  constraint 
since  we  do  not  have  enough  statistical  information  at  this  time.  Therefore,  an 
approximate  sound  speed  perturbation  covariance  matrix  is  generally  used.  When 
the  covariance  is  inexact,  the  estimator  is  suboptimal.  In  Chapter  IV,  we  vary  the 
assumed  correlation  length  (which  is  used  to  construct  the  sound  speed 
perturbation  covariance  matrix  for  the  estimator)  to  determine  the  sensitivity  of 
the  system  to  uncertainty  in  correlation  length.  We  also  study  the  effects  of 
experimental  random  noise  and  failure  of  array  elements  on  the  estimate  on  the 
estimate. 

In  Chapter  V,  we  present  a  summary  of  our  research,  along  with  results  and 
conclusions.  Furthermore,  we  propose  a  criterion  for  designing  estimators  for 
the  analysis  of  the  GSP  data.  Finally,  we  make  recommendations  for  improving 
our  research. 


II.  ACOUSTIC  FORWARD  MODELING 

A.  TRAVEL  TIME 

In  acoustic  tomography,  the  sound  speed  and  flow  fields  are  reconstructed 
from  travel  time  measurements.  The  corresponding  forward  problem  is, 
therefore,  to  find  eigenrays,  i.e.,  those  rays  which  emitted  by  the  source  that  are 
intercepted  by  the  receiver,  and  to  establish  the  relation  between  the 
measurement  and  the  unknown  fields.  Each  eigenray  has  a  unique  launch  angle, 
and  a  unique  path  through  the  ocean,  thereby  sampling  the  sound  speed  field  and 
flow  field,  at  different  locations  (Cornuelle,  1983,  pp.  41). 

The  geometric  approximate  travel  time  along  a  ray  path,  T(t),  can  be 
calculated  by  integrating  the  ray  slowness  along  the  path.  In  the  presence  of  a 
current   v(r,t),  T(t)  is  given  by 


I  ds 

JD„.u  c(r,t)  +  v(r,t) .  t 


T(t)  = 

Path   ~v"1'    '    'v"'v'      *  (2.1) 

where  c(r,t)  is  the  sound  speed  at  time  t  and  position  r,  v(r,t)  is  the  ocean 
current  velocity,  x  is  a  unit  vector  tangent  to  the  ray,  and  s  is  the  arc  length 
along  the  path.  The  travel  time  is  changed  by  the  sound  speed  perturbation  field, 
5c(r,t),  which  is  the  deviation  from  the  reference  sound  speed,  c0(r).  The 
reference  sound  speed  can  be  the  overall  space-time  average  sound  speed.  Thus 
the  travel  times  in  a  reciprocal  transmission  can  be  expressed  as 


Tf=T0  +  8Tf= 


and 


T    =  T„  +  8T  = 


f   ds  .     . 

/         c0(r,t)  +  8c(r,t)  +  v(r,t)  •  x 

"  path 

J        ds 

/         c0(r,t)  +  6c(r,t)-v(r,t)-x 

*  path 


(2.2) 


(2.3) 


where    the  superscripts  f  and  b  refer  to  forward  and  backward  transmissions, 
respectively,  and  8T  is  the  perturbation  of  the  reference  travel  time  T0. 

In  most  ocean  environments,  5c/c0  is  on  the  order  of  10    .  Thus  we  can 
approximately  linearize  the  reciprocal  travel  time  perturbations  as 


8T  =  - 


8c(r,t)  +  v(r,t)  •  x 


±ds 


co(0 


path 


and 


8T  =  - 


8c(r,t)  -  v(r,t)  •  x 


ds. 


c„(r) 


path 


Taking  the  sum  of  Eq.  (2.4)  and  Eq.  (2.5),  one  obtains: 


(2.4) 


(2.5) 


5t+=8tW 


8c(r,t) 
Co(r) 


ds, 


path 


while  taking  the  difference  of  Eqs.  (2.4)  and  (2.5),  one  obtains: 


(2.6) 


5T-=8Tf-5Tb  =  _     I        v^^ 
2  *  cfr) 

Path  (2.7) 

where  8T  ,  half  of  the  summation  of  the  forward  and  backward  travel  time 
perturbations,  is  linearly  related  to  the  sound  speed,  and  5T~,  half  of  the 
difference  of  the  forward  and  backward  travel  time  perturbations,  is  linearly 
related  to  the  current.  In  this  thesis  our  focus  is  on  the  estimation  of  the  sound 
speed  perturbations  only,  and  we  won't  be  dealing  with  Eq.  (2.7)  at  all. 

In  order  to  express  the  travel  time  in  a  vector  form  the  continuous  integral 
in  Eq.  (2.6)  was  discretized  by  dividing  the  ocean  into  small  boxes,  in  which  the 
sound  speed  perturbation  was  assumed  constant.  Figure  2-1  shows  the  horizontal 
transceiver  array  configuration  and  the  corresponding  box  geometry  in  a 
horizontal  slice. 

In  this  project  we  discretized  the  ocean  volume  by  500  boxes  (10x10  squares 
horizontally  and  5  layers  vertically);  the  limitation  of  computer  memory  space 
in  the  micro  VAX  dictated  this  decision. 

After  the  discretization  of  the  sound  speed  perturbation  field  into  a  vector 
5c.  the  vector  £T  containing  all  the  eigenray  travel  times  can  be  expressed  in  a 
matrix-vector  form  as 

£J  =  A£c,  (2.8) 

where 

A,.= =-    I    -rds, 


lJ    a  5cy  c; 


(2.9) 
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is  the  element  in  the  ith  row  and  y'th  column  of  the  matrix  A  and  is  equal  to  the 
integral  of  "  Vc2  along  a  segment  of  the  /th  acoustic  ray  path  in    the  yth  ocean 

box. 
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Figure  2.1:  The  Horizontal  Transceiver  Array  Configuration  and  the  Box 
Geometry  System   Used. 


B     ACOUSTIC  RAYS  IN    AN  INHOMOGENEOUS  MOVING 
MEDIUM 

A  numerical  ray-trace  algorithm  for  the  acoustic  rays  in  an  inhomogeneous 

moving  medium  was  used  in  our  study.  The  procedure  was  developed  by  Chiu 

(1985).  A  summary  of  the  theory  will  be  given  in  the  following  discussion. 


1 .     Ray  Path 

It  is  well  known,  in  acoustic  ray  theory,  that  the  acoustic  rays  in  a 
motionless  medium  with  space  dependent  index  of  refraction,  n(r),  is  described 
by  the  eikonal  equation: 

I  I2  2 

|VO(r)|  =n(r)  ,  (2.10) 

where  O(r)  is  the  eikonal  function  (i.e.,  acoustic  phase)  defining  the  wave  fronts, 
r  =  (x,y,z)  is  the  position  vector,  and  VO(r)  points  to  the  direction  of 
propagation  of  acoustic  field  energy.  However,  when  the  inhomogeneous 
medium  is  moving  with  a  velocity  given  by  v,  the  governing  eikonal  equation  for 
acoustic  wave  fronts  becomes  (Ugincius,  1970) 


2 

i2  2        ~    <S>(r) 

VO(r)    =n(r)  (1-v— — )  . 


c 


o 


(2.11) 


Ugincius  (1970)  has  derived  a  second-order  vector  differential  equation 
governing  the  ray  paths  based  on  the  eikonal  equation  Eq.  (2.11).  The  equation 
for  ray  paths  is 


4-  (Nr')  -  (rM«V')V  +  r'x(VxV)  =  VN, 

dsv  (2.12) 


where 


and 


N=|     ,     V=pv, 


=  V  l-z;2  +  (r'«v)2, 


n(S  -  r'*v) 
P  = 


(2.13) 


(2.14) 


d-1s' 


(2.15) 
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and  where  N  and  V  are  functions  of  both  the  position  r  of  a  point  as  well  as  the 
ray  direction  r'  through  that  point,  and  v  =  v/c  is  nondimensional  medium 
velocity  and  it  has  a  magnitude  equal  to  v.  Note  that,  we  have  used  '  and  "  to 
denote  the  first  and  second  derivatives  with  respect  to  s,  respectively. 

In  the  ocean,  the  current  velocity  is  small  compared  to  the  speed  of 
sound.  Therefore,  Eq.  (2.12)  can  be  simplified  by  eliminating  the  second  order 
terms  involving  v  in  N  and  V  so  that 

N  =  n     ,     V=  n(l  -r'»v)v  ,  (2.16) 

and 

(r".V)V-0.  (2>17) 

By  doing  so  Eq.  (2.12)  reduces  to 


j-  (nr')  +  r'x  Vxn(l  -  r'«v)v  ]  =  Vn. 


(2.18) 


By  replacing  r'  with  dr/ds  in  Eq.  (2.18)  we  have 


d  ,  dr.      dr  |_     ...     dr      . 
ds^(nds-)+ds-xVxn(1-ds-v)v 


=  Vn. 


(2.19) 

In  this  thesis  we  neglected  the  ray  curvature  in  the  horizontal  plane  and 
let  the  reference  sound  speed  profile  to  be  a  function  of  depth  z  only.  This 
restriction  implies  that  the  reference  ray  paths  are  confined  to  lie  on  vertical 
slices  normal  to  the  xy  plane  (Ziomek,  1985,  pp.  232).  In  Figure  2-2  a  ray  path 
with  such  a  restriction  is  shown.  On  a  vertical  slice  (i.e.,  a  Range-Depth  Rz 
plane),  the  equations  governing  the  planar  rays  on  this  plane  are 
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and 


d  ,  dR,      d ,      . 
ds      ds      ds      ^ 


d  ,  dz,      d ,     . 


tf  = — R 

dR 


(2.20) 


z  =  —  z, 
9z 


ds     ds      ds 
where  R  and  z  are  the  unit  vector  in  R  and  z  direction,  respectively. 


(2.21) 


z  (depth) 

Figure  2-2:  The  Ray  Path  Restricted  in  a  Vertical  Slice  Normal  to  xy 
Plane. 

Since  the  reference  sound  speed  profile  is  taken  to  be  a  function  of 

depth  only,  we  can  simplify  the  differential  equation  (Eq.  (2.20))  to  get 


dR     cos80-P0+n(z)i;R 


ds 


n(z) 


=  cos(0) 


(2.22) 


or 


dz 
dR 


-V 


n(z) 


cos8-i70+n(z)rR 


-  1   =  tan(8). 


(2.23) 

where  vQ  is  the  medium  speed  at  position  (R0,  z0)  in  the  ray  path  direction.  As 
shown  in  Figure  2-2  80  and  8  are  the  initial  launch  angle  at  position  (R0,  z0)  and 
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ray  angle  at  position  (R,  z),  respectively;  Eq.  (2.23)  can  be  integrated  in  either 
the  R  or  z  directions  to  get  the  ray  paths.  In  our  raytracing  we  divided  each 
range  connecting  each  of  the  transceiver  pairs  in  the  Greenland  Sea  into  1,000 
steps  and  integrated  Eq.  (2.23)  in  the  R  direction  for  the  corresponding  depth 
values.  Integrating  in  range  gives,  for  each  ray  path,  depth  of  the  trajectory  as  a 
function  of  range.  The  function  has  a  one-to-one  correspondence. 

2.     Numerical  4th  Order  Runge-Kutta  Integration  Method 

An  accurate  ray  path  can  be  calculated  using  a  well  known  numerical 
4th  order  Runge-Kutta  integration  method.  From  Eq.  (2.23).  the  integral  for 
which  we  need  to  compute  is: 


z=z0± 


4 


n(z) 


-  1    dR 


cos0o-  v0+  n(z)i'R 

(2.24) 

The    numerical  4th  order  Runge-Kutta  integration  method  is  given  by 

the  following  formulae  (Gerald,  1989,  pp.  358): 


J9     Rn+h  1 

Zn+1  =  Zn  +  Jr         f  [R.  Z(R)1  dR  =  Zn  +  £  (k  !+  2k2  +  2k3  +  k4)  , 


where 


and 


(2.25.1) 


f[R,z(R)]  =  ±/v/ 

n(z) 

2 
-1     , 

cos0o-D  0+n(z)i;  R 

(2.25.2) 

kx  =  h  f(Rn  ,  zj, 

(2.25.3) 

k2=hf(Rn+ih,zn  +  ik1), 

(2.25.4) 

k3  =  hf(Rn  +  ^h 

,  Zj,  +  2  k2), 

(2.25.5) 
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k4=hf(Rn+h,  zn+k3).  (2.25.6) 

A  step  size  h  has  been  selected  to  limit  the  numerical  error  to  a 

tolerable  size.  The  global  error  (i.e.,  discretization  error)  accumulated  along  the 

entire  interval  R  by  the  4th  order  Runge-Kutta    method  is  (Gerald,  1989,  pp. 

358) 

o(h5)  =  f^f4M]     ,      o<^<h.  (226) 

The  error  cannot  be  exactly  determined  because  the  position  (R,  z)  =  [^  ,  z(  Q]  is 
an  unknown  with  §  bounded  by  the  interval  [0,  h].  A  standard  way  to  determine 
whether  the  z  values  are  sufficiently  accurate  is  to  compare  the  value  computed 
using  a  step  size  of  h  with  the  value  calculated  using  the  half  of  h.  If  this  gives 
only  a  change  of  negligible  magnitude,  the  results  are  accepted;  if  not,  the  step  is 
halved  again  until  the  results  are  satisfactory  (Gerald,  1989,  pp.  358). 
3.     Turning  Point 

As  an  acoustic  ray  travels  through  an  ocean  volume  it  will  be  refracted 
upward  or  downward.  When  the  ray  angle  goes  to  zero,  the  ray  will  start  to 
bend  up  or  down.  The  position  of  this  zero  ray  angle  point  is  called  the  turning 
point. 

When  a  ray  is  traced  to  a  position  which  is  less  than  a  step  size  away 
from  a  turning  point,  we  use  a  linear  gradient  approximation  to  calculate  the 
depth  of  the  ray  at  the  terminal  of  that  step.  This  approximation  near  a  turning 
point  is  needed  because  the  function  for  which  we  integrate  will  no  longer  be  the 
same  after  the  turning  point.  There  will  be  a  sign  change  in  the  function.  The 
curvature  K  can  be  used  to  calculate  the  ray  path  around  turning  point  in  the 
linear  gradient  case.  Ugincius  (1970)  has  derived  the  following  equation  of 
curvature: 
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1 

K  =  — 
C 


dc  ,      dR,     dv  ,~dR       1N 

d^(l7-dF)+dF(2d^-1) 


(2.27) 

Let  the  gradient  of  sound  speed   dc/dz  be  gc  and  the  gradient  of  medium  speed 
dv/dz  be  gv  Eq.  (2-27)  becomes 


K  = 


1 


gc  (v  -  cosG)  +  g v  (2cos0  v  -  1) 


(2.28) 


Since  gv  and  v    are  very  small,  the  term  involving  the    product  of  gj, 
and  v  is  negligible.  Thus  Eq.  (2.28)  can  be  approximated  by 


K  = 


glv  -cosQ)-gv) 


The  radius  of  curvature  ^,is  the  reciprocal  of  k  ,  i.e., 


(2.29) 


(2.30) 


£=r y 

gjfr  -cosQ)-gv) 

Under  the  assumption  of  constant  gradients  near  a  turning  point,  the 
local  radius  !^.is  also  a  constant,  which  result  in  a  circular  ray  path,  locally.  If 
the  radius  ^,  were  negative,  the  radius  would  curve  upward  and  vice  versa 
(Kinsler,  1982,  pp.  401-402).  In  the  following  discussion  ^,will  be  referred  to  as 
the  magnitude  of  the  radius  of  curvature. 

The  geometries  of  segments  of  ray  paths  through  turning  points  are 
shown  in  Figure  2-3  and  Figure  2-4,  where  Qx  is  the  ray  angle  before  the  turning 
point  and  62  is  the  ray  angle  after  the  turning  point.  The  corresponding  depth 
increment  can  be  equated  as 


z/+l-zJ=±(^cos02"  ^COsOj), 


(2.31) 
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step  size  h    — ►] 


R 


z  (depth) 
Figure  2-3:  The  Geometry  of  a  Downward  Turning  Ray  Path. 


-    R 


z  (depth) 
Figure  2-4:  The  Geometry  of  an  Upward  Turning  Ray  Path. 
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where 

n T 

^cosQ2  =  y  ^K,-  (h-^sine,)  ,  (2  32) 

and  "+"  is  for  downward  curve  while  "-"  is  for  upward  curve.  From  Eq.  (2.31) 
the  depth  of  the  turning  point  is  simply  given  by 

zTO=z/±^(l-cos91).  (2  33) 

4.    Surface   Reflection 

Due  to  the  upward  refracting  nature  of  the  Greenland  Sea  sound 
channel,  surface  reflection  of  rays  is  of  importance.  The  ray  path  is 
discontinuous  at  the  point  of  reflection  and  the  equation  for  rays  Eq.  (2.24)  is  not 
valid  here.  We  need  to  develop  a  procedure  to  calculate  the  local  path  trajectory 
near  a  reflection  point  in  an  other  way.  The  same  linear  gradient  assumption,  as 
used  for  the  calculation  near  a  turning  point  is  applicable  here  again. 

Similar  to  Eq.  (2.31),  the  depth  increment  in  a  range  step  within  which 
a  surface  reflection  occurs  is 

zJ+1-zJ  =  ±(!^cos02-  ^cose^.  p  34^ 

In  the  case  of  convex  path  segment  (figure  2-5),  ^,cos02  is  given  by 


^.cos92  =  V  1L  -  (2^sin9  -  ^,sin9 ,  -  h  ) 

and 

!£sin9  =  y  ^  -  (!£cos9  {  +  Z}   . 

In  the  case  of  a  concave  path  segment  (Figure  2-6),  ^.cos02  is  given  by 


(2.35) 


(2.36) 
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Figure  2-5:  The  Geometry  of  a  Surface  Reflected  Convex  Path  Segment. 


Figure  2-6:  The  Geometry  of  a  Surface  Reflected  Concave  Path  Segment. 
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^COS0. 


=Y^ 


(2^sin0-  ^.sin9j  +  h  ) 


(2.37) 


and 


=v7^ 


^sin0  =  'Y  9L-  (^cosej-Z^)   .  (2  38) 

Having  all  the  equations  coded  in  FORTRAN,  predictions  of  ray  paths  as  well  as 
eigenray  finding  were  possible. 
5.     Eigenrays   Finding 

Since  the  sound-speed  profile  in  the  central  Greenland  Sea  gyre  is  very 
nearly  adiabatic  below  the  surface  layers,  a  linear  and  range-independent  sound 
speed  profile  was  selected  to  represent  the  reference  state  c0(z),  as  shown  in 
Figure  2-7  (left).  We  then  use  this  reference  state  to  calculate  the  ray  pattern  of 
transmission  in  the  acoustic  forward  problem. 


RAY  TRACE  BETWEEN  ARRAY  «1  &  »2 


1  .45     1  .50  0 
C  (Ka/s) 


40  60 

RANGE  I  Km) 


Figure  2-7:  Typical    Sound  Speed  Profile  and  Ray  Path  with  The  Ray  Angle 
from  -15°  to  0°  in  The  Greenland  Sea  Project  Area. 
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In  Figure  2-7  (right)  the  ray  paths  shown  are  the  typical  propagation 
pattern  in  the  arctic  environment.  The  rays  with  fewer  loops,  that  reach  the 
deeper  layers,  are  the  faster  paths  for  the  acoustic  energy. 

In  order  to  search  for  rays  that  reach  the  receiver  (i.e.,  eigenrays)  we 
need  to  shoot  rays  with  a  range  of  launch  angles.  Given  an  angular  interval 
within  which  all  possible  eigenrays  lie,  we  shall  be  able  to  determine  the 
eigenrays  by  looking  for  the  intersections  between  the  arrival  depth  curve  and 
the  receiver  depth  line  as  illustrated  in  Figure  2-8.  Three  eigenray  patterns 
associated  with  3  different  ranges  are  shown  in  Figure  2-9.  We  only  choose  6  to 
7  out  of  as  many  as  30  eigenrays  that  sample  the  ocean  from  the  surface  to  1  km 
depth.  Figure  2-10  shows  the  selected  eigenrays  for  the  three  ranges.  There  are  a 
total  of  91  eigenrays  will  be  used  to  conduct  our  sensitivity  study. 
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Figure  2-8:  Arrival  Depth  Curve  for  Eigenrays  Finding. 
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Figure  2-9:  Typical  Eigenray  Pattern  for  Three  Different  Range  as  in  (a) 
Range  =  123.6  km,  (b)  Range  =  105.2  km,  and  (c)  Range  =  200.1  km. 
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Figure  2-10:  Typical     Selected  Eigenray  Pattern  for  Three  Different    Range 
as  in  (a)  Range  =  123.6  km,  (b)  Range  =  105.2  km,  and  (c) 
Range  =  200.1  km. 
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III.     STOCHASTIC  INVERSE  METHOD 

In  the  forward  problem  the  ray  travel  times  are  modeled  using  Eq.  (2.9).  In 
this  chapter,  we  discuss  the  inversion  of  Eq.  (2.9)  using  a  Gauss-Markoff 
estimator  to  obtain  an  optimal  estimate  of  the  three  dimensional  sound  speed  field 
in  the  ocean  volume  monitored  by  the  acoustic  array.  The  estimation  error  and 
system  resolution  will  be  analyzed  in  the  next  chapter  in  an  effort  to  quantify  the 
performance  of  the  array. 

A.   ESTIMATION  OF  SOUND  SPEED  FIELDS. 

1.      The  Gauss-Markoff  Estimator 

In  the  Gauss-Markoff  stochastic  inverse  method,  both  the  data  £T  and 
the  unknown  field  £c  are  assumed  to  be  random  vectors.  The  forward  model 
relating  the  data  and  the  unknown  field  is 

§T  =  A£c  +  e,  (3  1) 

where  e  is  the  experimental  noise  which  corrupts  the  travel  time  measurement;  e 
is  assumed  to  be  uncorrected  with  £T  and  §c. 

The  system  Eq.  (3.1)  is  highly  underdetermined,  and  thus  without 
additional  constraints  other  than  just  the  data,  the  system  admits  infinite  number 
of  solutions  for  £c_.  In  stochastic  inverse  methods  a  unique  optimal  linear 
estimate  of  the  unknown  parameters,  £c,  can  be  constructed  by  incorporating  a 
priori  knowledge  of  the  parameters  in  a  covariance  matrix. 

The  well  known  Gauss-Markoff  estimator  is  chosen  here  because  it 
requires  no  knowledge  of  probability  densities.  This  so  called  distribution  free 
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property  is  the  most  important  feature  of  the  Gauss- Markoff  estimator  (Liebelt, 
1967,  pp.  136).  The  estimate  £c_  satisfies  a  minimum  mean  square  error 
criterion 

(  I  £c  -  £c  I    )  =  minimum.  n  9^ 

Following  Liebelt  (1967)  and  Chiu  (1987),  the  Gauss-Markoff  estimate 
is  given  by 

*  T    — 1 

§c.  =  C£A  Ce  £T,  £3.3) 

where 


c£= 


-1 

:^+atc:1a 


(3.4) 
is  the  covariance  matrix  of  the  error  £  =  £c-  £c  in  the  estimate,  and  Cg  and  C^ 

are  the  covariance  matrices  of  noise  e  and  the  unknown  parameters  §c, 
respectively.  The  construction  of  the  estimate  requires  finding  the  inverse  of  the 

-1  T         -1 

matrix  C^  +   A   C£A. 

The  trading  of  system  resolution  for  stability  in  the  optimal  estimate  can 

i       i 

be  revealed  by  a  singular  value  decomposition  of  the  matrix  C  2AC  2   such  that 


C  2AC2R 

e  5c 


T 

=  U  A  V  , 

(3.5) 


where  the  diagonal  elements  A,,  of  the  matrix  A  are  the  associated  singular  values, 
and  the  columns  u,  and  v,  of  U  and  V  are  the  left  and  right  singular  vectors, 
respectively. 

The    matrix,    Eq.    (3.5),    is    the    operator    associated    with    a 
nondimensionalized  version  of  the  forward  model: 
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C  2£1  = 


C7ACL 

e  £c 


c&& +  CS> 


(3.6) 


which  is  transformed  from  Eq.  (3.1).  Using  Eq.  (3.5),  Chiu  et  al  (1987)  have 
shown  that  the  minimum  mean  square  error  estimate  can  be  expressed  as 


c^a=v 


I  +  A 


-l 


or  equivalently  as 


ck  =  i 

«•  =  i 


X: 


&f 


X:  +  1 


U     C  2£T 

e 


11,-C/filUi 


(3.7) 


(3.8) 


The  vectors  v,  are  the  base  vectors  occupying  the  solution  domain  such 
that  any  solution  can  be  expressed  by  a  weighted  sum  of  these  vectors.  A  singular 
value  A,,  much  small  than  one  is  associated  with  a  singular  vector  y^  that  models 
a  highly  unstable  and  oscillatory  function.  The  linear  estimator  downweights 
these  oscillatory  function  to  stabilize  the  estimate. 

B.   ERROR  AND  RESOLUTION 

In  the  previous  section,  we  introduced  an  estimator  which  solves  the  inverse 
problem.  In  this  section,  we  derive  measures  which  are  used  to  quantify  system 
performance. 

1.     Error  of  the  Estimate 

The  error  covariance  matrix  of  the  estimate  was  expressed  in  Eq.  (3.4), 
which  can  also  be  equated,  using  Eq.  (3.5),  as 


C£  "C5c 


I-VAI  +  A      AV 


-l 


'&' 


(3.9) 


The  diagonal  terms  of  this  matrix  are  the  mean  square  errors  of  the   estimate  at 
each  box. 


.A 


The  error  of  the  estimate,  §c-£c,  has  two  components,  bias  and 


.A 


random  error.  The  bias,  b  =  (§c)  -  §c,  results  from  an  insufficient  number  of 

data,  and  the  random  error,  A(£c  )  =    5c  -  (&c >  is  caused  by    the  random 
experimental  noise.  Since  these  two  components  are  statistically  independent,  C£ 

can  be  expressed  as 


where 


or  equivalently, 


C£=  (bb    )  +  CA(5£> 


T     -  1 

^A(3c)=C£A    Cg  AC£ 


(3.10) 


(3.11) 


cA(52)  -  C^ 


,-1 

2  2 


-1 


V  1+A      A    I  +  A        V 


'5c' 


(3.12) 


An  expression  for  (bb    )  can  be  obtained  by  letting  the  random  error 
covariance  CA(Ss) approaches  zero,  i.e.  by  letting  all  the  eigenvalues  ^approach 

infinity  in  Eq.  (3.9).  The  result  is 


<<>bT>=c; 


£c  I 


I- vv 


(3.13) 


2.     Resolution 

We  define  a  symmetric,  nxn  matrix 


(i  +  A2)    i 


R  =  VA\I  +  A  /    AV     .  (3  14) 

as  the  resolution  matrix.  Using  the  definition  of  R,  Eq.  (3.9)  becomes 
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c£=cL[i-R]  c 


5c  l  j      sc  •  (3  15) 

If  the  resolution  matrix  is  an  identity  matrix  (that  is,  R=I),  then  the 
error  is  zero.  On  the  other  hand,  if  the  resolution  matrix  is  not  an  identity 
matrix,  for  example,  if  it  has  nonzero  diagonal  elements  then  the  error  is 
nonzero.  This  implies  that  the  mean  square  error  is  intimately  related  to  system 
resolution.  The  ith  row  of  the  resolution  matrix,  {rj  ,  is  defined  as  the 
resolution  kernel  of  the  ith  box  and  it  describes  how  much  neighboring  boxes 
contribute  to  the  error  of  the  estimate  in  the  ith  box,  or  how  well  the  field  at  the 
ith  box  can  be  resolved  (Menke,  1984,  pp.  61-68).  Figure  3-1  illustrates  the 
structure  of  a  typical  resolution  matrix. 


Figure  3-1:  Plot  of  Selected  Rows  and  Columns  of  The  Resolution  Matrix 
R. 


If  the  ith  resolution  kernel  has  a  single  spike  centered  at  the  diagonal, 
then  the  ith  box  is  perfectly  resolved.  If  the  peak  is  very  broad,  the  ith  box  is 
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poorly  resolved  and  has  a  large  error  due  to  the  fact  that  the  estimate  at  the  ith 
box  is  an  average  of  the  neighboring  field.  In  an  other  word  the  energy  in  the 
estimate  is  spread  to  neighboring  boxes  due  to  the  poor  resolving  power  of  the 
system. 

Two  measures  of  the  resolution  of  the  estimator  have  been  proposed. 
One  is  the  so  called  resolution  "spread"  (Miller,  1989,  pp.  333),  which  measures 
the  difference  between  the  resolution  matrix  and  an  identity  matrix  via  the 
expression 

n        n  2 

spread  =  £    X  (Ry- ly)  > 

i  =  lj  =  l  (3.16) 

The  second  measure  is  the  so  called  minimum  resolution  length  (Chiu, 

1987),  which  measures  the  local  resolution  at  each  box  location.  The  resolution 

length  is  defined  to  be  the  distance  at  which  the  resolution  energy  falls  to  half  of 

its  peak  value.  To  be  more  precise  the  minimum  resolution  length  at  the  ith  box 

is  defined  as  the  square  root  of  the  second  central  moment  of  energy  distribution 

in  the  ith  resolution  kernel  ({r,}  ).   The  minimum  resolution  lengths  in  the  three 

spatial  directions  may  be  expressed  as 


2r?(jx,jyjz) 


Xx(ix,iy7iz)  =  /i  /  £  £  £  [  (jx  -  ix)  dx  ] 


jx    jy     jz  « 


(3.17.a) 


/  nx  iw    nz  2  r -(jxjyjz) 

^(ix,iy,iz)  =  /i  /  £  £  £  [  (jy  -  iy)  dy  ]   J^L^Ll , 


J*  ^   Jz  '  (3.17.b) 

and 
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a*i    ■    m         /v^Vrr       -u    ,2  rfQx  jy,iz) 
1\ix,iy,iz)  =  s\l  Z,Z,2*[(jz-iz)dz]   — 

Y    i*  jy  iz  ' 

where  E,is  the  total  energy  of  the  /th  resolution  kernel  r/. 


(3.17.C) 


(3.17.d) 


j»  jy   jz 

In  Eq.  (3.17.a-b)  ^  and  .7^  are  the  minimum  horizontal  resolution  lengths  in  the 
x  and  y  directions,  respectively,  while  V  is  the  minimum  vertical  resolution 
length.  Physically,  the  minimum  resolution  lengths  determine  the  minimum  eddy 
size  that  can  be  resolved  adequately  by  the  monitoring  system.  Note  that  we  have 
expanded  the  row  index  i  and  column  index  j  into  three  dimensional  box  indices, 
(ix,  iy,  iz)  and  (jx,  jy,  jz),  in  Eqs.  (3.17.a-d).  The  expansion  was  done  according 
to  the  following  equations. 

i  =  (ix-l)xnyxnz  +  (iy-l)xnz  +  iz  =  1,  2 n,  (3. 18. a) 

j  =  (jx-l)xnyxnz  +  (jy-l)xnz  +  jz  =  1,  2 n  (3.18.b) 

with 

n  =  nxxnyxnz,  (3. 1 8.c) 

where  ix,  jx  =  1,  2,  nx,  is  box  indices  in  the  x  direction;  iy,  jy  =  1,2, 

ny,  is  box  indices  in  the  y  direction;  iz,  jz  =  1,2, nz,  is  box  indices  in  the 

z  direction.  With  nx  =  10,  ny  =  10,  and  nz  =  5,  the  total  number  of  boxes  is  n  = 
500. 

From  the  above  discussion  in  this  chapter,  It  has  been  found  that  the 
RMS  error  and  resolution  length  do  not  depend  on  the  experimental  data  §T,  but 
only  on   the  sound  speed  perturbation  covariance  matrix  C$c,  travel  time  error 
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covariance  matrix  C^,  and  the  transfer  function  A  (which  depends  on  the  array 
geometry  and  the  number  of  rays). 
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IV.     RESULTS  OF  SENSITIVITY  STUDY 

In  the  last  chapter,  we  outlined  two  ways  to  evaluate  the  performance  of  the 
Greenland  Sea  tomography  array  by  examining  the  statistics  of  the  error  in  the 
estimate  and  the  resolution  of  the  system.  As  mentioned  the  RMS  error  of  the 
estimate  and  the  system  resolution  measures  depend  on  the  covariance  matrix  of 
sound  speed  perturbation  C^,  the  covariance  matrix  of  noise  C£,  the  array 

geometry,  the  number  of  eigenrays  but  not  the  data  5T.  Therefore,  even  without 
using  actual  or  synthetic  measurements,  the  performance  of  the  array  and  ray 
path  geometry  can  be  evaluated  given  the  covariances.  By  varying  the  ocean 
correlation  length  (which  determines  C^),  the  rms  value  of  noise  (which 

determines  Ce),  and  the  array  geometry  (which  determines  A),  we  have 
examined  the  changes  in  system  performance  as  these  ocean  and  acoustic 
parameter  vary.  The  results  are  discussed  in  Section  IV-B. 

By  incorporating  statistical  information  concerning  the  covariance  of  the 
field,  the  indeterminacy  of  the  unknown  field,  5c,  is  eliminated.  Since  the 
covariance  matrix  C^  is  the  a  priori     information  that  we  supply  to  the 

estimator,  we  are  particularly  interested  in  determining  the  sensitivity  of  the 
estimate  to  the  uncertainties  in  the  correlation  lengths  of  the  sound  speed 
perturbation  field.  This  sensitivity  study  was  accomplished  through  inversions  of 
synthetic  data  generated  in  the  computer.  The  results  are  discussed  in  Section  IV- 
C. 

In  the  following  section,  we  first  discuss  the  method  we  use  to  simulate  sound 
speed  perturbation  fields  and  travel  time  data  in  the  computer. 
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A.   COMPUTER  SIMULATION  OF  MESOSCALE  SOUND  SPEED 
FIELDS  AND  TRAVEL  TIME  DATA 

The  temperature,  sigma-t,  and  sound  speed  profiles  obtain  by  a  CTD 
(Conductivity,  Temperature,  and  Density)  cast  in  the  Greenland  Sea  (Worcester 
and  Howe,  1989)  is  shown  in  Figure  4-1.  We  have  superimposed  on  this  data  a 
linear  profile  (dash  line)  which  we  have  chosen  as  the  reference  sound  speed 
profile.  In  the  simulation  work,  we  take  that  the  perturbations  of  sound  speed 
occur  only  in  the  water  column  shallower  than  1,000  m.  This  is  not  a  bad 
assumption  as  indicated  by  CTD  profile  obtained  by  Worcester  and  Howe.  For 
simplicity,  we  also  neglect  ocean  currents  in  our  analysis  and  work  with  forward 
transmission  paths  only. 

In  the  reconstruction  of  the  sound  speed  perturbation  field,  the  correlation 
covariance  function  of  the  field  (or  its  discretized  version,  i.e.,  the  covariance 
matrix  C^.)  needs  to  be  specified.  We  assume  that  the  sound  speed  perturbation 

field  is  homogeneous  and  has  a  gaussian  shape  so  that  we  can  specify  the 
correlation  between  the  field  at  two  different  points  in  an  analytical  form  as 


2 

Cov(Ax,Ay,Az)  =  o*§^e 


\2    /     \2    /     -2 

Ax     +   A>M  +   Az 
L.  U  L, 


(3.8) 

where  Ax  is  the  horizontal  separation  between  the  two  points  in  the  x-direction, 
Ay  is  the  horizontal  separation  in  the  y-direction,  and  Az  is  the  vertical 
separation.  The  correlation  length,  Lx,  Ly,  and  Lz,  determine  the  correlation 
scale  of  the  field.  TABLE  4-1  summarizes  the  various  sets  of  correlation  lengths 
(Lx,  Ly,  and  Lz)  that  we  used  to  simulate  the  sound  speed  perturbation  fields  of 
different  scales.  For  all  the  simulated  fields,  we  use  <38c  -  5  m/s. 
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Figure  4-1:  The  CTD  Data  of  Mooring  #1  From  the  Deployment  Cruise 
(Worcester  and   Howe,   1989). 
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TABLE  4-1 

:  THE  SI 

Ocean  volume 

Lx  (km) 

Eddy204 

20 

Eddy304 

30 

Eddy404 

40 

THE  SIMULATED  OCEAN  VOLUME. 


Ly  (km) 
20 
30 
40 


Lz  (km) 
0.4 
0.4 
0.4 


In  order  to  perform  simulation  inversions  in  the  computer,  we  need  to 
generate  a  set  of  travel  time  data  for  input  to  the  estimator.  A  normal  random 
deviate  generator  is  used  to  generate  realizations  of  a  random  process  having  a 
specified  covariance.  This  random  deviate  generator  is  used  to  simulate  the  sound 
speed  perturbation  and  noise  fields.  The  simulated  fields  are  then  combined  using 
the  model  (given  by  Eq.  (2.9))  to  give  the  simulated  travel  time  data.  A  block 
diagram  of  the  process  is  shown  in  Figure  4-2. 


■6c 


Random 
deviator 
generator 

Random 
deviator 
generator 

-5T 


Figure  4-2  :  The  Block  Diagram  of  Travel  Time  Generation. 

The  generated  travel  time  data  with  additive  random  noise  are  the  input  to 
the  estimator.  An  estimate  of  the  sound  speed  field,  £c  is  the  inverse  result, 
which  depend  on  the  matrices  C^,  C£  and  A  in  addition  to  the  data  themselves. 

The  sensitivity  of  the  system  can  thus  be  evaluated  for   changes  in  correlation 
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length  of  the  sound  speed  perturbation  field,  changes  in  noise  level,  changes  in 
array  geometry  due  to  mooring  failure,  etc.  To  simulate  noise-free 
measurements,  we  use  a  very  small  noise  level  with  an  rms  value  of  ae  =  0.1  ms. 
The  optimal  estimate  is  obtained  when  the  correlation  length  used  for 
inversion  is  exactly  the  same  as  that  actually  present  in  the  ocean  volume  (that  is, 
when  the  a  priori  covariance  C^  is  correct).  In  fact,  because  the  a  priori 

covariance  is  never  exactly  known,  optimal  estimates  are  generally  hard  to 
obtain.  The  quality  of  suboptimal  estimates  can  be  evaluated  by  studying  the 
effect  of  correlation  length  uncertainties  ALX,  ALy,  and  ALZ  on  the  performance 
of  the  estimator. 

B.    SYSTEM  PERFORMANCE 
1.     RMS  Error  Analysis 

The  local  RMS  error  as  a  function  of  box  index  j  in  the  estimate  (that 

is,  the  RMS  error  at  thejth  box)  can  be  obtained  by  evaluating  the  square  root  of 
the  v'th  element  along  the  diagonal  of  the  estimate  error  covariance  matrix  C£. 

This  local  RMS  error  gives  a  picture  of  how  the  errors  are  distributed  spatially. 
The  square  root  of  the  spatial  average  of  the  mean  square  errors  at  each  of  the 
boxes,  aE,  was  calculated  to  get  an  idea  of  how  much  the  local  error  varies  over 
the  whole  ocean  volume  being  studied. 

Figure  4-3  shows  the  contour  plots  of  the  local  RMS  estimate  error  for 
one  set  of  system  parameters  in  each  of  the  five  vertical  layers.  The  RMS  error 
for  this  case  is  approximate  from  1  m/s  to  2.5  m/s  (or  20%  to  50%  compared  to 
the  5  m/s  signal  level)   and  shows  gradual  dependence  on  horizontal  or  vertical 
position    inside  the  perimeter  of  the  array.    The  error  outside   the   array's 
perimeter  increases  rapidly  as  distance  from  the  array  increases.  Obviously,  we 
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Figure  4-3:  RMS  Error  (m/s)  at  Each  Layer.  The  contour  level  is  0.5  m/s. 
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cannot  obtain  accurate  estimate  of  the  sound  speed  perturbations  for  points 
outside  the  array.  There  is  little  difference  between  the  RMS  errors  in  the 
individual  layers  because  each  vertical  layer  are  well  sampled  by  the  eigenrays 
whose  turning  points  occupy  every  layer.  Because  of  the  rather  weak  vertical 
dependence,  in  the  following  discussion  we  present  only  the  analysis  for  the  first 
layer. 

The  local  RMS  error  maps  for  the  first  layer  for  different  horizontal 
correlation  lengths  are  shown  in  Figure  4-4.  The  system  parameters  are  identical 
to  those  of  Figure  4-3,  except  that  the  horizontal  correlation  length  is  varied 
from  20  km  to  60  km  in  a  10  km  step.  We  see  that  a  wider  co variance  matrix 
(i.e.,  a  longer  correlation  length)  results  in  a  lower  estimate  error.  This 
observation  is  not  surprising,  since  a  longer  correlation  length  means  that  there 
is  more  gradual  variation  in  sound  speed  perturbation  with  position,  and  that 
neighboring  boxes  are  more  correlated.  An  increased  correlation  reduces  the 
number  of  degrees  of  freedom  in  the  solution,  and  thus  giving  a  better 
determined  solution. 

In  Figure  4-5,  the  maps  are  generated  using  the  system  parameters 
identical  to  those  used  to  generate  Figure  4-3,  except  that  the  vertical  correlation 
length  is  varied  from  0.2  km  to  0.6  km  in  a  0.1  km  step  at  a  fixed  horizontal 
correlation  length  of  40  km.  There  is  less  change  in  RMS  error  over  the  entire 
range  of  correlation  lengths  for  this  case  than  for  the  case  shown  in  Figure  4-4. 
This  is  because  the  eigenrays  sample  the  vertical  layer  more  adequately  than  the 
horizontal  sections. 

The  effect  of  an  increased  (or  decreased)  noise  level  on  the  estimate  are 
shown  in  Figure  4-6.  A  higher  noise  level  results  only  in  a  slightly  higher  RMS 
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Figure  4-4:   Estimate  RMS  Error  (m/s).  The  contour  level   is  0.5  m/s.  The 
vertical  correlation  length  Lz  is  fixed  at  0.4  km. 
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Figure  4-5:   Estimate  RMS  Error  (m/s).  The  contour  level   is  0.5  m/s.  The 
horizontal  correlation  length  Lx  and  Ly  is  fixed  at  40  km. 
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Figure  4-6:  Estimate  RMS  Error  at  Different  Noise  Levels.  The  contour 
level   is  0.5  m/s. 
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error.  The  estimate  is  thus  not  as  strongly  affected  by  changes  in  noise  as  it  is  by 
changes  in  the  ocean  correlation  length. 

Figure  4-7  shows  the  effect  of  array  element  failure  on  the  RMS  error 
(the  array  configuration  is  superimposed  on  the  contour  plots  to  give  an  idea  of 
how  the  estimate  depends  on  the  array  configuration).  Within  the  perimeter  of 
the  remaining  "good"  elements,  the  RMS  errors  are  only  about  25%  (0.5  m/s) 
higher  than  those  given  by  the  full  array  (2.0  m/s).  However,  the  errors  outside 
the  area  covered  by  the  good  elements  increase  rapidly  to  100%.  The  failure  of 
array  elements  has  a  very  pronounced  effect  on  the  estimate,  especially  when  two 
or  more  elements  fail.  The  system  essentially  loses  its  ability  to  give  accurate 
estimates  in  areas  that  do  not  have  rays  passing  through  them. 

A  single  measure  of  RMS  error  can  be  evaluated  by  calculating  a  spatial 
"average"  of  the  individual  local  errors.  We  use  oE  for  such  a  global  measure 

and  it  is  computed  by  assuming  the  individual  mean-square  errors  at  all  the  box 

locations  and  then  taking  the  square  root  of  that  sum.  Figures  4-8  and  4-9  show 
oEas  a  function  of  the  ocean  horizontal  and  vertical  correlation  lengths, 
respectively,  for  various  array  geometries  and  noise  levels.  In  all  cases,  GE  was 
found  to  decrease  with  increasing  Lxand  Ly  (see  Figure  4-8).  However,  aEis  not 
as  sensitive  to  changes  in  Lz  as  it  is  to  changes  in  Lx  and  Ly  (see  Figure  4-9).  This 
is  due  to  the  fact  that  the  vertical  structure  is  sampled  adequately  whereas  the 
horizontal  structure  is  not.  In  view  of  this,  in  the  following  discussion,  we  will 
only  discuss  the  sensitivity  of  the  results  of  our  resolution  analysis  with  the 
vertical  correlation  length  fixed  at  0.4  km. 
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Figure   4-7:    Estimate   RMS    Error   (m/s)    in    Mooring   Failure   Cases.   The 
contour  level  is  0.5  m/s.  Lx=  Ly  =  40  km,  and  Lz  =  0.4  km. 
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Figure  4-8:  Square  Root  of  Spatial  Average  Mean  Square  Error  vs. 
Horizontal  Correlation  Length. 
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Figure  4-9:  Square  Root  of  Spatial  Average  Mean  Square  Error  vs.  Vertical 
Correlation  Length. 
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2.     Resolution  Analysis 

Two  measures  of  system  resolution  were  defined  in  the  last  chapter. 
One  measure  is  the  minimum  resolution  length,  which  gives  the  local  resolution 
at  each  box.  The  minimum  resolution  length  is  essentially  the  size  of  the  smallest 
ocean  feature  which  can  be  resolved  by  the  an-ay.  A  large  minimum  resolution 
length  indicates  a  poor  resolving  power.  Figure  4-10  shows  the  minimum 
resolution  length  ^  in  the  x  direction  as  a  function  of  Lx  and  Ly.  We  see  that 
for  Lx  and  Ly  greater  than  30  km,  !tfx  is  relatively  constant  over  the  interior  of 
the  array  and  in  most  of  the  region  of  interest.  On  the  other  hand,  if  Lx  and  Ly 
are  less  than  or  equal  to  30  km,  ^  varies  heavily  on  both  x  and  y  and  becomes 
very  large  in  regions  containing  no  y-oriented  ray  paths.  The  behavior  of  9L,  the 
minimum  resolution  length  in  the  y  direction,  is  analogous  to  the  behavior  of  ^ 
and  is  shown  in  Figure  4-11.  In  general,  the  average  minimum  resolution  length 
is  approximately  30  km  inside  the  monitored  region  as  shown  in  Figure  4-10  and 
Figure  4-11. 

Figure  4-12  shows  the  behavior  of  Mx  and  9L  as  noise  level  change 

while  the  horizontal  and  vertical  correlation  lengths  are  fixed  at  40  km  and  0.4 
km,  respectively.  We  see  that  varying  the  noise  level  has  little  effect  on  ^  and 
iHy.  Figure  4-13  shows  how  the  failure  of  various  array  elements  affects  2/"x.  A 
diagram  of  the  array  configuration  is  superimposed  on  each  plot  in  Figure  4-13 
to  show  the  connection  between  array  geometry  and  'H^.  Despite  the  element 
failures,  ^  remains  relatively  constant  inside  the  area  covered  by  the  remaining 
array  elements,  but  it  increases  rapidly  as  we  move  outward  from  the  perimeter. 
Although  there  is  some  resolving  power  outside  the  array  perimeter,  it  is  highly 
inadequate. 
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Figure  4-10:  Horizontal  Minimum  Resolution  Length  0{x.  The  contour  level 
is  10  km.  Lz  is  fixed  at  0.4  km. 
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Figure  4-11:  Horizontal  Minimum  Resolution  Length  9{y.  The  contour  level 


is  10  km.  Lz  is  fixed  at  0.4  km. 
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Figure  4-12:  Horizontal  Minimum  Resolution  Lengths  at  different 
levels.  The  contour  level  is  10  km.  Lx  =  Ly  =  40  km,  and  Lz 
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The  second  measure  of  system  resolution  is  the  resolution  spread  (Eq. 
3.16),  which  is  the  square  of  the  Frobenius  norm  of  the  difference  between  the 
resolution  matrix  and  the  identity  matrix.  The  resolution  spread  is  thus  a  single 
scalar  quantity  which  describes  the  resolution  of  an  estimator  over  an  entire 
region.  Figure  4-14  shows  resolution  spread  as  a  function  of  horizontal 
correlation  length  for  each  of  several  noise  levels  and  array  geometries.  We  see 
that  the  resolution  spread  does  not  significantly  depend  on  changes  in  correlation 
length,  except  when  noise  level  is  increased.  That  is,  system  resolution  is  a  lot 
more  sensitive  to  changes  in  correlation  length  when  the  system  is  not  noise  free. 
Generally,  the  resolution  spread  depends  most  strongly  on  the  array  geometry  at 
any  fixed  correlation  length. 
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Figure  4-14:  Resolution  Spread  vs.  Horizontal  Correlation  Length. 
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C.    SENSITIVITY  OF  INVERSE  SOLUTION 

Synthetic  travel  time  perturbation  data  £X  (refer  to  Figure  4-2  for  the 
generation  of  £T)  were  used  to  reconstruct  the  sound  speed  perturbation  fields  £c 
simulated  in  the  computer.  The  sensitivity  of  the  estimate  to  uncertainty  in  the 
correlation  length  specified  for  inversion  is  examined  in  this  section  by 
comparing  the  suboptimal  estimates  to  the  optimal  estimates.  The  optimal 
estimates  are  those  derived  using  an  exact  covariance  of  the  field  whereas  the 
suboptimal  ones  are  consequences  of  inexact  covariances. 

Figure  4-15  shows  the  simulated  field  Eddy404  (refer  to  TABLE  4-1  for  the 
simulation  parameters),  its  optimal  estimate,  the  associated  RMS  error  as  well  as 
the  difference  between  the  estimate  and  the  simulated  field.  The  estimate  error 
§c-£c  is  low  (from  0  to  ±2  m/s,  which  is  from  0%  to  40%  compared  to  a  signal 
level  of  5  m/s)  over  most  of  the  area  inside  the  array.  The  error  is  larger  in  the 
left  edge  and  corners.  Figures  4-16  compares  the  estimate  errors  of  some 
suboptimal  estimates  generated  using  correlation  lengths  differ  from  the  true 
one.  The  difference  between  the  assumed  horizontal  correlation  lengths  for 
inversion  and  the  true  lengths  are  denoted  by  ALX  and  ALy  in  the  figure. 
Obviously,  the  effect  of  a  positive  correlation  length  uncertainty  seems  less 
harmful  than  that  of  a  negative  one  when  the  actual  correlation  length  is  40  km. 

Figure  4-17  shows  the  effect  of  uncertainty  in  the  noise  level  (i.e.,  in  the 
noise  covariance  matrix  CJ  on  the  estimate.  The  maps  on  the  left  are  the 
estimates  of  the  Eddy404  field  with  various  noise  level  uncertainties  A<7g.  On  the 
right  the  associated  errors  are  displayed.  Generally,  a  higher  noise  level 
uncertainty  gives  a  higher  estimate  error,  although  the  differences  in  the  errors 
are  quite  minimal. 
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Figure  4-15:  A  simulated  field,  its  optimal  estimate,  estimate  RMS  error, 
and  difference  between  the  optimal  estimate  and  the  simulated  field.  Units 

are  in  m/s.  Lx=  Ly  =  40  km,  Lz  =  0.4  km,  and  ALX=  ALy  =  ALZ=  0. 
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Figure    4-16:    Simulated    Edd>404    field,    suboptimal    estimates,    and    the 
corresponding  true  error  fields.    The  contour  level  is  4  m/s  and  2  m/s  for 

£c.  and  §£  -  §£,  respectively.  Lx  =  Ly  =  40  km,  and  Lz  =  0.4  km. 
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Figure  4-16:  (Continued  from  last  page). 
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Figure  4-17:   Estimation   (m/s)  with   Noise  Level   Uncertainty.  The  contour 

level  is  4  m/s  and  2  m/s  for  ££  and  $£  -  &£,  respectively.  Lx  =  Ly  = 
40  km,  and  Lz  =  0.4  km. 
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To  gain  insight  into  the  global  performance  of  the  estimator,  we  calculated 
G^,  the  square  root  of  the  spatial  average  of  the  local  squared  errors.  Figure  4-18 
shows  o^  as  a  function  of  horizontal  correlation  uncertainty  for  each  of  the  three 
simulated  ocean  volumes  having  horizontal  correlation  length  of  20,  30  and  40 
km,  respectively. 


(m/s) 


d     eddy202-h 


Figure  4-18:  Square  Root  of  Spatial  Average  Square  Error  vs. 
Horizontal  Correlation  Length  Uncertainty. 


It  is  interesting  to  note  that,  for  assumed  horizontal  correlation  lengths  of  Lx 
and  Ly  both  equal  to  1  km*,  the  error  is  3.4  m/s  for  all  the  three  cases  estimated. 
Furthermore,  a^  is  less  sensitive  to  a  positive  correlation  length  uncertainty  than 


*  Correlation  length  of  zero  implies  no  a  priori  information.  Our  software 
does  not  permit  the  use  of  zero  correlation  length,  so  we  used  correlation  length 
of  1  km  to  approximate  the  case  of  no  a  priori  information. 
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to  a  negative  correlation  length  uncertainty  when  the  actual  ocean  correlation 
length  is  long  (e.g.  longer  than  20  km,  see  the  curves  associated  with  ocean  fields 
Eddy304  and  Eddy404),  while  it  is  more  sensitive  to  a  positive  uncertainty  than 
to  a  negative  uncertainty  when  the  actual  correlation  length  is  short  (e.g.  20  km, 
see  the  curve  associate  with  Eddy 204). 

In  an  ocean  with  small  eddies,  Figure  4-18  suggests  that  using  a  correlation 
length  longer  than  the  actual  one  could  result  in  a  large  increase  of  estimate 
error.  Whether  or  not  the  estimate  is  sensitive  to  an  inexact  and  longer 
correlation  length  used  really  depends  on  the  resolution  of  the  system.  We  see 
from  Figure  4-10  that  the  minimum  horizontal  resolution  length  of  the 
Greenland  Sea  array  is  approximately  30  km,  almost  regardless  of  what  the 
correlation  length  of  the  field  is.  This  implies  that  the  array  is  unable  to  resolve 
ocean  features  smaller  than  30  km.  In  other  words,  an  ocean  field  with  a 
correlation  length  shorter  than  30  km  (e.g.  20  km)  can  not  be  monitored 
adequately.  In  that  case  any  additional  a  priori  information,  even  though  wrong, 
is  totally  absorbed  by  the  estimator  for  solution  construction,  thus  leading  to  an 
uncertainty  sensitive  estimate.  In  contrast,  for  an  ocean  field  with  a  correlation 
length  longer  than  30  km,  the  array  turns  into  an  adequate  system.  The 
corresponding  estimator  now  has  the  ability  to  reject  extraneous  information. 
The  result  is  an  uncertainty  insensitive  estimate. 

TABLE  4-2  shows  the  value  of  a^  for  each  suboptimal  estimate,  as  well  as 
the  percent  difference  in  a^  between  each  suboptimal  estimate  and  the  optimal 

estimate  for  the  ocean  field  Eddy404.  By  selecting  a  value  for  the  maximum 
acceptable  percent  difference  in  a^  between  the  suboptimal  estimators  and  the 

optimal  estimator,  we  can  compute  from  the  curves  in  Figure  4-18  the  maximum 
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allowable  uncertainty  in  correlation  length.  TABLE  4-3  shows,  for  each 
simulated  field,  the  allowable  range  of  uncertainty  in  correlation  length  if  the 
difference  between  a^  for  the  suboptimal  and  optimal  estimators  is  not  to  exceed 

ten  percent. 

TABLE  4-2  :  THE  ESTIMATE  ERROR  AND     %  DIFFERENCE  TO 
THE  OPTIMAL  ESTIMATE  FOR  EDDY404. 


ALX,  ALy 

(km) 

°S 

% 

difference 

-30 

3.30 

75.53 

-20 

2.50 

32.98 

-10 

2.06 

9.57 

10 

1.88 

0.00 

20 

1.99 

9.63 

30 

2.10 

11.70 

TABLE  4-3  :  TEN   %  DIFFERENCE  ALLOWABLE  UNCERTAINTY 
RANGE. 

Lx,  Ly  (km)  Allowable  ALX  and  ALy  Range  (km) 

20  km  -20.0-15.6 

30  km  -11.8-23.7 

40  km  -10.4-25.8 

TABLE  4-3  quantitatively  confirms  our  earlier  observation  from  Figure  4- 
18  that  the  estimate  error  is  less  sensitive  to  a  positive  uncertainty  in  correlation 
length  than  to  a  negative  one  in  ocean  volumes  containing  large  structures,  and 
the  result  is  reversed  for  an  ocean  volume  containing  small  structures. 
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V.  CONCLUSION 

A.  ESTIMATOR  PERFORMANCE 

In  order  to  obtain  a  unique  estimate  of  the  sound  speed  perturbation  field,  a 
priori  information  in  the  form  of  a  sound  speed  perturbation  covariance  matrix 
is  used  in  the  Gauss-Markoff  estimator.  As  discussed  in  Chapter  4,  given  that  the 
sound  speed  perturbation  field  has  a  gaussian  shape  correlation,  the  optimal 
estimate  for  the  field  is  obtained  when  the  assumed  correlation  length  (i.e.,  the 
correlation  length  used  to  calculate  the  input  covariance  matrix  C^)  is  equal  to 

the  actual  correlation  length  present  in  the  ocean.  A  primary  goal  of  this  thesis 
is  to  evaluate  the  effect  on  our  estimator  when  the  assumed  correlation  length  is 
not  equal  to  the  actual  one.  As  this  happens,  the  estimate  becomes  suboptimal. 
When  the  actual  correlation  length  is  not  exactly  known,  it  has  been  suggested  by 
Cornuelle  (1985)  and  Chiu  (1987)  to  use  a  conservative  assumption  {i.e.,  a  small 
correlation  length)  so  as  not  to  "assume  too  much"  about  the  sound  speed 
perturbation  field. 

From  our  simulation  study  we  found  that  the  optimal  estimate  for  the  sound 
speed  perturbation  field  typically  has  an  RMS  error  between  1  and  2  m/s  (i.e., 
20%  to  40%  comparing  to  the  5  m/s  signal  level).  For  the  suboptimal  case,  the 
estimate  is  actually  less  sensitive  to  a  positive  correlation  length  uncertainty  {i.e., 
the  assumed  correlation  length  is  longer  than  the  actual  one)  than  to  a  negative 
uncertainty  when  the  actual  correlation  length  is  longer  than  30  km.  On  the 
other  hand,  when  the  actual  correlation  length  is  shorter  than  30  km  (for 
example  20  km)  the  estimate  becomes  more  sensitive  to  a  positive  uncertainty 
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than  to  a  negative  uncertainty.  The  reason  for  this  switching  of  behavior  at  a 
correlation  length  of  30  km  has  something  to  do  with  the  resolution  of  the 
Greenland  Sea  array.  In  our  resolution  analysis  we  found  that  the  Greenland  Sea 
array  has  a  resolution  length  of  about  30  km.  Therefore,  for  an  ocean  field  with 
a  correlation  length  small  than  30  km  (for  example  20  km),  the  use  of  a 
correlation  length  longer  than  the  actual  one  for  inversion  would  ingest 
extraneous  information  into  an  information  hungry  estimator.  This  estimator 
basically  accepts  all  the  wrong  information.  It  is  thus  preferable  to  be 
conservative  and  use  a  correlation  length  which  is  likely  to  have  a  negative 
uncertainty  when  the  system  resolution  is  inadequate  for  measuring  the  expected 
scale.  On  the  other  hands,  if  we  know  from  resolution  analysis  that  the  resolution 
is  adequate  for  a  particular  ocean  region,  a  positive  correlation  length 
uncertainty  is  acceptable  in  this  case. 

By  specifying  the  maximum  acceptable  percent  difference  in  o^  (which  is 

essentially  a  spatial  average  of  local  estimate  errors)  between  the  suboptimal  and 
optimal  estimators,  we  arrive  at  one  possible  criterion  for  designing  the 
estimator.  If  this  difference,  in  the  case  of  an  expected  ocean  correlation  length 
of  about  40  km,  is  required  to  be  less  than  ten  percent  (for  example),  the 
estimator  can  tolerate  any  correlation  length  uncertainty  between  -10.4  and  25.8 
km,  a  spread  of  35.6  km. 

We  have  also  shown  that  a  higher  noise  level  results  in  basically  unchanged 
RMS  error  and  resolution.  This  result  suggests  that  the  estimate  error  is  not 
dominated  by  the  random  error,  but  rather  by  the  bias  error  arising  from  the 
fact  that  the  tomography  problem  is  underdetermined.  This  insensitivity  to 
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random  noise  is  one  benefit  of  using  a  Gauss-Markoff  estimator,  which  always 
tries  to  minimize  the  effect  of  random  noise  (Chiu  et  aly  1987). 

The  failure  of  array  elements  has  a  very  pronounced  effect  on  the  RMS 
error  and  resolution,  especially  if  two  elements  fail.  The  comparison  of  cases 
simulating  the  failure  of  one  or  more  elements  indicates  that  the  RMS  error 
becomes  very  large  and  resolution  becomes  very  poor  in  areas  that  no  longer 
have  rays  passing  through  them.  However,  in  regions  still  containing  acoustic 
rays,  the  RMS  error  is  only  25%  higher  than  that  of  the  full  array  case. 

B.    RECOMMENDATIONS  FOR  FUTURE  IMPROVEMENT 

All  numerical  simulations  were  performed  on  a  DEC  Micro  VAX.  Each  run 
required  roughly  five  and  one-half  hours  of  CPU  time.  A  machine  of  greater 
computational  power  and  larger  memory  size  would  allow  us  to  divide  the  ocean 
volume  into  finer  spatial  meshes  for  improved  analysis. 

We  can  derive  statistical  information  concerning  the  vertical  ocean  structure 
from  historical  data  using  the  empirical  orthogonal  function  (EOF)  approach  of 
Cornuelle  (1983,  pp.  139).  EOF  analysis  has  been  widely  applied  in  research 
fields  other  than  tomography.  In  EOF  analysis,  the  vertical  structure  of  the  ocean 
is  represented  by  a  set  of  orthogonal  vectors.  These  vectors  can  be  derived  from 
a  singular  value  decomposition  of  a  matrix  containing  historical  data  from 
hydrographic  surveys.  The  EOF  method  gives  a  priori  information  about  the 
vertical  structure  which  is  more  realistic  than  that  provided  by  the  Gaussian 
shape  correlation  used  in  this  study,  and  is  recommended  for  use  when  the  actual 
tomographic  data  from  the  Greenland  Sea  Project  become  available. 
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